
/** 计算儒略日 */
export function JulDay (time: Date): number {
    const ut = time.getUTCHours() + time.getUTCMinutes()/60 + time.getUTCSeconds()/3600;
    let y = time.getUTCFullYear();
    let m = time.getUTCMonth() + 1; // 月份从0开始，所以需要加1
    const d = time.getUTCDate();
    if (m<=2) {m=m+12; y=y-1}
    return Math.floor(365.25*(y+4716)) + Math.floor(30.6001*(m+1)) + d - 13 -1524.5 + ut/24.0;
}

function frac(X: number): number {
    X = X - Math.floor(X);
    if (X<0) X = X + 1.0;
    return X;
}

/** 十进制小时 */
type DecHour = number;

/** 十进制度 */
type DecDegree = number;

/** 弧度 */
type Radian = number;

export function HoursMinutesSeconds(time: DecHour): string {
    // 把十进制小时转换为"HH:MM:SS"格式
    const hours = Math.floor(time);
    const minutes = Math.floor((time - hours) * 60);
    const seconds = Math.floor(((time - hours) * 60 - minutes) * 60);

    return `${String(hours).padStart(2, '0')}:${String(minutes).padStart(2, '0')}:${String(seconds).padStart(2, '0')}`;
}

function GM_Sidereal_Time (jd: number): number {
	var t_eph, ut, MJD0, MJD;
	MJD = jd - 2400000.5;
    MJD0 = Math.floor(MJD);
	ut = (MJD - MJD0)*24.0;
	t_eph  = (MJD0-51544.5)/36525.0;
	return  6.697374558 + 1.0027379093*ut + (8640184.812866 + (0.093104 - 0.0000062*t_eph)*t_eph)*t_eph/3600.0;
}

/**
 * 计算本地恒星时
 * @param jd 儒略日
 * @param longitude 经度
 * @returns
 */
export function LM_Sidereal_Time (jd: number, longitude: DecDegree): DecHour {
	var GMST = GM_Sidereal_Time(jd);
	var LMST =  24.0*frac((GMST + longitude/15.0)/24.0);
  	return LMST;
}

// ========== 行星位置计算 ==========

export const deg2rad = (deg: DecDegree): Radian => deg * Math.PI / 180;
export const rad2deg = (rad: Radian): DecDegree => rad * 180 / Math.PI;
export const rad2hour = (rad: Radian): DecHour => rad * 12 / Math.PI; // 弧度转为十进制小时

export const planetsConst = [
    {name: "水", n: 4.092317, M0: 174.795, a: 0.38710, e: 0.20563, i: 7.005, omega: 29.125, Omega: 48.331, A: 0.37073, color: "#d0d0d0"},
    {name: "金",   n: 1.602136, M0:  50.416, a: 0.72333, e: 0.00677, i: 3.395, omega: 54.884, Omega: 76.680, A: 0.72330, color: "#f0d0a0"},
    {name: "地",   n: 0.985608, M0: 357.529, a: 1.00000, e: 0.01671, i: 0.000, omega: 288.064, Omega: 174.873, A: 0.99972, color: "#00a0f0"}, // 地球
    {name: "火",    n: 0.524039, M0:  19.373, a: 1.52368, e: 0.09340, i: 1.850, omega: 286.502, Omega: 49.558, A: 1.51039, color: "#ff6040"},
    {name: "木", n: 0.083056, M0:  20.020, a: 5.20260, e: 0.04849, i: 1.303, omega: 273.867, Omega: 100.464, A: 5.19037, color: "#f8d870"},
    {name: "土",  n: 0.033371, M0: 317.021, a: 9.55491, e: 0.05551, i: 2.489, omega: 339.391, Omega: 113.666, A: 9.52547, color: "#f0c080"},
    {name: "天",  n: 0.011698, M0: 141.050, a: 19.21845, e: 0.04630, i: 0.773, omega: 98.999, Omega: 74.006, A: 19.17725, color: "#a0e0f0"},
    {name: "海", n: 0.005965, M0: 256.225, a: 30.11039, e: 0.00899, i: 1.770, omega: 276.340, Omega: 131.784, A: 30.10796, color: "#5080ff"},
    {name: "冥",   n: 0.003964, M0:  14.882, a: 39.543, e: 0.2490, i: 17.140, omega: 113.768, Omega: 110.307, A: 37.09129, color: "#a07070"},
]

/** 地球黄赤交角 */
const earthE: Radian = deg2rad(23.439292);

export function getPlanetPosition(time: Date) {
    const jd = JulDay(time);
    const days = jd - 2451545.0; // J2000 epoch

    // 计算日心坐标系下的行星位置
    const heliocentricPositions = planetsConst.map(planet => {
        let M = planet.M0 + planet.n * days; // Mean anomaly
        M = M % 360; // Normalize to [0, 360)
        M = M * Math.PI / 180; // Convert to radians
        const niu = M + planet.e; // True anomaly
        const r = planet.A / (1 + planet.e * Math.cos(niu)); // Distance in AU

        const Omega = planet.Omega * Math.PI / 180; // Longitude of ascending node
        const omega = planet.omega * Math.PI / 180; // Argument of periapsis
        const i = planet.i * Math.PI / 180; // Inclination
        const x = r * (Math.cos(Omega) * Math.cos(niu + omega) - Math.sin(Omega) * Math.sin(niu + omega) * Math.cos(i));
        const y = r * (Math.sin(Omega) * Math.cos(niu + omega) + Math.cos(Omega) * Math.sin(niu + omega) * Math.cos(i));
        const z = r * (Math.sin(niu + omega) * Math.sin(i));
        const distance = Math.sqrt(x * x + y * y + z * z); // Distance in AU
        const elon = Math.atan2(y, x) // Ecliptic longitude
        const elat = Math.asin(z / distance); // Ecliptic latitude
        return {x, y, z, distance, elon, elat, name: planet.name};
    });

    // 计算地心坐标系下的行星位置
    const geocentricPositions = heliocentricPositions.map((pos, index) => {
        // 地心直角坐标系
        if (index === 2) { // Earth
            return {
                ra: 0,
                dec: 0,
            }
        }
        const x = pos.x - heliocentricPositions[2].x; 
        const y = pos.y - heliocentricPositions[2].y;
        const z = pos.z - heliocentricPositions[2].z;
        const distance = Math.sqrt(x * x + y * y + z * z); // 地心距离
        const elon = Math.atan2(y, x); // 地心黄经
        const elat = Math.asin(z / distance); // 地心黄纬

        // 计算赤经赤纬
        const ra = Math.atan2(Math.sin(elon) * Math.cos(earthE) - Math.tan(elat) * Math.sin(earthE), Math.cos(elon));
        const dec = Math.asin(Math.sin(elat) * Math.cos(earthE) + Math.cos(elat) * Math.sin(earthE) * Math.sin(elon));
        return {ra, dec,}
    });

    // 整合返回
    return geocentricPositions.map((_, index) => {
        const hp = heliocentricPositions[index];
        const ep = geocentricPositions[index];
        return {
            name: hp.name,
            hDistance: hp.distance, // 日心距离
            hElon: hp.elon, // 日心黄经
            hElat: hp.elat, // 日心黄纬
            hx: hp.x, // 日心直角坐标系X
            hy: hp.y, // 日心直角坐标系Y
            hz: hp.z, // 日心直角坐标系Z
            ra: ep.ra, // 赤经
            dec: ep.dec, // 赤纬
        };
    });
}

/**
 * 将赤经赤纬转换为方位角和高度
 * @param siderealTime 本地恒星时，单位十进制小时
 * @param lat 纬度，单位弧度
 * @param ra 赤经，单位弧度
 * @param dec 赤纬，单位弧度
 * @returns { azi: 方位角, alt: 高度 }
 */
export function re2azi(siderealTime: DecHour, lat: Radian, ra: Radian, dec: Radian){
    const ha = (siderealTime - rad2hour(ra)) * 15 * Math.PI / 180; // 时角，单位弧度
    const z = Math.acos(Math.sin(lat) * Math.sin(dec) + Math.cos(lat) * Math.cos(dec) * Math.cos(ha));
    const azi = Math.atan2(Math.sin(ha), Math.cos(ha) * Math.sin(lat) - Math.tan(dec) * Math.cos(lat)); // 方位角，单位弧度
    const alt = Math.PI / 2 - z; // 高度，单位弧度
    return { azi, alt };
}